{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Segmentation: Thresholding and Edge Detection</h1>\n",
    "\n",
    "In this notebook our goal is to estimate the location and radius of spherical markers visible in a Cone-Beam CT volume.\n",
    "\n",
    "We will use two approaches:\n",
    "1. Segment the fiducial using a thresholding approach, derive the sphere's radius from the segmentation. This approach is solely based on SimpleITK.\n",
    "2. Localize the fiducial's edges using the Canny edge detector and then fit a sphere to these edges using a least squares approach. This approach is a combination of SimpleITK and scipy/numpy.\n",
    "\n",
    "Note that all of the operations, filtering and computations, are natively in 3D. This is the \"magic\" of ITK and SimpleITK at work.\n",
    "\n",
    "The practical need for localizing spherical fiducials in CBCT images and additional algorithmic details are described in:\n",
    "Z. Yaniv, \"Localizing spherical fiducials in C-arm based cone-beam CT\", *Med. Phys.*, Vol. 36(11), pp. 4957-4966."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(SimpleITK)\n",
    "\n",
    "source(\"downloaddata.R\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the volume and look at the image (visualization requires window-leveling)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "spherical_fiducials_image <- ReadImage(fetch_data(\"spherical_fiducials.mha\"))\n",
    "Show(spherical_fiducials_image, \"spheres\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After looking at the image you should have identified two spheres. Now select a Region Of Interest (ROI) around the sphere which you want to analyze. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "roi1 = list(c(280,320), c(65,90), c(8, 30))\n",
    "roi2 = list(c(200,240), c(65,100), c(15, 40))\n",
    "mask_value = 255\n",
    "\n",
    "# Select the ROI\n",
    "roi = roi1\n",
    "# Update the R ROI, SimpleITK indexes are zero based, R indexes start at one  \n",
    "r_roi = lapply(roi, function(x) x+1)\n",
    "    \n",
    "# Create the mask image from an R array    \n",
    "amask <- array(0, spherical_fiducials_image$GetSize())\n",
    "xs <- r_roi[[1]][1]:r_roi[[1]][2]\n",
    "ys <- r_roi[[2]][1]:r_roi[[2]][2]\n",
    "zs <- r_roi[[3]][1]:r_roi[[3]][2]\n",
    "amask[xs, ys, zs] <- mask_value\n",
    "\n",
    "mask <- Cast(as.image(amask), \"sitkUInt8\")\n",
    "mask$CopyInformation(spherical_fiducials_image)\n",
    "Show(LabelOverlay(Cast(IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767, \n",
    "                                          windowMaximum=-29611), \n",
    "                       \"sitkUInt8\"), \n",
    "                   mask, opacity=0.5))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholding based approach\n",
    "\n",
    "Our region of interest is expected to have a bimodal intensity distribution with high intensities belonging to the spherical marker and low ones to the background. We can thus use Otsu's method for threshold selection to segment the sphere and estimate its radius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set pixels that are in [min_intensity,otsu_threshold] to inside_value, values above otsu_threshold are\n",
    "# set to outside_value. The sphere's have higher intensity values than the background, so they are outside.\n",
    "\n",
    "inside_value <- 0\n",
    "outside_value <- 255\n",
    "number_of_histogram_bins <- 100\n",
    "mask_output <- TRUE\n",
    "\n",
    "labeled_result <- OtsuThreshold(spherical_fiducials_image, mask, inside_value, outside_value, \n",
    "                                number_of_histogram_bins, mask_output, mask_value)\n",
    "\n",
    "# Estimate the sphere radius and location from the segmented image using the LabelShapeStatisticsImageFilter.\n",
    "label_shape_analysis <- LabelShapeStatisticsImageFilter()\n",
    "label_shape_analysis$SetBackgroundValue(inside_value)\n",
    "label_shape_analysis$Execute(labeled_result)\n",
    "cat(\"The sphere's center is: \", label_shape_analysis$GetCentroid(outside_value), \"\\n\")\n",
    "cat(\"The sphere's radius is: \",label_shape_analysis$GetEquivalentSphericalRadius(outside_value),\"mm\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Edge detection based approach\n",
    "\n",
    "In this approach we will localize the sphere's edges in 3D using SimpleITK. We then compute a least squares sphere that optimally fits the 3D points using scipy/numpy. The mathematical formulation we use is as follows:\n",
    "\n",
    "Given $m$ points in $\\mathbb{R}^n$, $m>n+1$, we want to fit them to a sphere such that\n",
    "the sum of the squared algebraic distances is minimized. The algebraic distance is:\n",
    "$$\n",
    "\\delta_i = \\mathbf{p_i}^T\\mathbf{p_i} - 2\\mathbf{p_i}^T\\mathbf{c} + \\mathbf{c}^T\\mathbf{c}-r^2\n",
    "$$\n",
    "\n",
    "The optimal sphere parameters are computed as:\n",
    "$$\n",
    "[\\mathbf{c^*},r^*] = argmin_{\\mathbf{c},r} \\Sigma _{i=1}^m \\delta _i ^2\n",
    "$$\n",
    "\n",
    "setting $k=\\mathbf{c}^T\\mathbf{c}-r^2$ we obtain the following linear equation system ($Ax=b$):\n",
    "$$\n",
    "\\left[\\begin{array}{cc}\n",
    "      -2\\mathbf{p_1}^T & 1\\\\\n",
    "      \\vdots        & \\vdots \\\\\n",
    "     -2\\mathbf{p_m}^T & 1\n",
    "     \\end{array}\n",
    "\\right]\n",
    "\\left[\\begin{array}{c}\n",
    "      \\mathbf{c}\\\\ k\n",
    "      \\end{array}\n",
    "\\right] =\n",
    "\\left[\\begin{array}{c}\n",
    "      -\\mathbf{p_1}^T\\mathbf{p_1}\\\\\n",
    "      \\vdots\\\\\n",
    "      -\\mathbf{p_m}^T\\mathbf{p_m}\n",
    "      \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "The solution of this equation system minimizes $\\Sigma _{i=1}^m \\delta _i ^2 = \\|Ax-b\\|^2$.\n",
    "\n",
    "Note that the equation system admits solutions where $k \\geq\n",
    "\\mathbf{c}^T\\mathbf{c}$. That is, we have a solution that does not\n",
    "represent a valid sphere, as $r^2<=0$. This situation can arise in\n",
    "the presence of outliers.\n",
    "\n",
    "Note that this is not the geometric distance which is what we really want to minimize and that we are assuming that there are no outliers. Both issues were addressed in the original work (\"Localizing spherical fiducials in C-arm based cone-beam CT\").\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create a cropped version of the original image.\n",
    "sub_image = spherical_fiducials_image[r_roi[[1]][1]:r_roi[[1]][2],\n",
    "                                      r_roi[[2]][1]:r_roi[[2]][2],\n",
    "                                      r_roi[[3]][1]:r_roi[[3]][2]]\n",
    "\n",
    "# Edge detection on the sub_image with appropriate thresholds and smoothing.\n",
    "edges <- CannyEdgeDetection(Cast(sub_image, \"sitkFloat32\"), \n",
    "                            lowerThreshold=0.0, \n",
    "                            upperThreshold=200.0, \n",
    "                            variance = c(5.0, 5.0, 5.0))\n",
    "\n",
    "# Get the 3D location of the edge points\n",
    "edge_indexes <- which(as.array(edges)==1.0, arr.ind=TRUE)\n",
    "# Always remember to modify indexes when shifting between native R operations and SimpleITK operations\n",
    "physical_points <- t(apply(edge_indexes - 1, MARGIN=1,\n",
    "                           sub_image$TransformIndexToPhysicalPoint))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually inspect the results of edge detection, just to make sure. Note that because SimpleITK is working in the\n",
    "physical world (not pixels, but mm) we can easily transfer the edges localized in the cropped image to the original."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "edge_label <- Image(spherical_fiducials_image$GetSize(), \"sitkUInt8\")\n",
    "edge_label$CopyInformation(spherical_fiducials_image)\n",
    "e_label <- 255\n",
    "dev_null <- apply(physical_points, \n",
    "                  MARGIN=1, \n",
    "                  function(x, img, label) img$SetPixel(img$TransformPhysicalPointToIndex(x),label), \n",
    "                  img=edge_label, \n",
    "                  label=e_label)\n",
    "    \n",
    "Show(LabelOverlay(Cast(IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767, windowMaximum=-29611),\n",
    "                       \"sitkUInt8\"), \n",
    "                  edge_label, opacity=0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup and solve linear equation system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A <- -2 * physical_points\n",
    "A <- cbind(A, 1)\n",
    "b <- -rowSums(physical_points^2)\n",
    "x <- solve(qr(A, LAPACK=TRUE), b)\n",
    "cat(\"The sphere's center is: \", x[1:3], \"\\n\")\n",
    "cat(\"The sphere's radius is: \", sqrt(x[1:3] %*% x[1:3] - x[4]), \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now, solve using R's linear model fitting. We also weigh the edge points based on the gradient magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gradient_magnitude = GradientMagnitude(sub_image)\n",
    "grad_weights = apply(edge_indexes-1, MARGIN=1, gradient_magnitude$GetPixel)\n",
    "\n",
    "df <- data.frame(Y=rowSums(physical_points^2), x=physical_points[, 1],\n",
    "                 y=physical_points[, 2], z=physical_points[, 3])\n",
    "fit <- lm(Y ~ x + y + z, data=df, weights=grad_weights)\n",
    "\n",
    "center <- coefficients(fit)[c(\"x\", \"y\", \"z\")] / 2\n",
    "radius <- sqrt(coefficients(fit)[\"(Intercept)\"] + sum(center^2))\n",
    "cat(\"The sphere's center is: \", center, \"\\n\")\n",
    "cat(\"The sphere's radius is: \", radius, \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## You've made it to the end of the notebook, you deserve to know the correct answer\n",
    "\n",
    "The sphere's radius is 3mm. With regard to sphere location, we don't have a the ground truth for that, so your estimate is as good as ours. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
